p <- 0.1

sample_sizes <- c(10,20,50,200)

par(mfrow=c(2,2))

for( n in sample_sizes) {
  x <- 0:n
  y <- dbinom(x, size=n, prob=p)
  
  plot(x, y, type="h", lwd=2, col="blue",
       main=paste("Binomial Distribution (n =", n, ", p =", p, ")"),
       xlab="Number of successes", ylab="Probability",
       ylim=c(0, max(y) * 1.1)
       )
  
  points(x, y, pch=19, col="red")
}

par(mfrow=c(1,1))

weights <- c(2,5,8,14,18)
proportion <- c(0.15, 0.35, 0.2, 0.15, 0.15)

m_result <- sample(weights, size=10, prob=proportion, replace=TRUE)

sm_func <- function(sz) {
  sample(weights, size=sz, prob=proportion, replace=TRUE)
}

mean(m_result)
sd(m_result)

barplot(table(m_result), 
        main="Sampled Weights Distribution",
        xlab="Weights", ylab="Frequency",
        col="lightblue", border="blue")


rep_result <- replicate(10000, mean(sm_func(10)))

m1 <- mean(rep_result)
sd1<- sd(rep_result)

hist(rep_result, breaks=seq(0, 16, by=0.5), 
     main="Histogram of Sampled Weights",
     xlab="Weights", ylab="Frequency",
     col="lightgreen", border="darkgreen", prob = TRUE)

# lines(density(rep_result), col="red", lwd=2)

curve(dnorm(x, mean=m1, sd=sd1), 
      from=min(rep_result), to=max(rep_result), 
      col="blue", lwd=2, add=TRUE)


